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Abstract: We present an image reconstruction algoritlim for tiie Inverse Conductivity Problem based on reformulating 
the problem in terms of integral equations. We use as data the values of injected electric currents and of the corresponding 
induced boundary potentials, as well as the boundary values of the electrical conductivity. 

We have used a priori information to find a regularized conductivity distribution by first solving a Fredholm integral equation 
of the second kind for the Laplacian of the potential, and then by solving a first order partial differential equation for 
the regularized conductivity itself. Many of the calculations involved in the method can be achieved analytically using the 
eigenfunctions of an integral operator dcBned in the paper. 
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1. Introduction 

It is well-known that the inverse conductivity problem is an extremely ill-posed, non-linear inverse 
problem, and there has been much interest in determining the class of conductivity distributions that 
can be recovered from the boundary data, as well as in the development of related reconstruction 
algorithms. The interest in this problem has been generated by both difficult theoretical challenges and 
by the important medical, geophysical and industrial applications of this problem. Much theoretical 
work has been related to the approach of Calderon concerning the bijection between the Neumann-to- 
Dirichlet operator, which relates the distribution of the injected currents to the boundary values of the 
induced electrical potential and the conductivity inside the region [1-4] . The reconstruction procedures 
that have been proposed include a wide range of iterative methods based on formulating the inverse 
problem as a nonlinear optimisation problem. These techniques are quite demanding computationally 
particularly when addressing the three dimensional problem. This concern has encouraged the search 
for reconstruction algorithms which reduce the computational demands either by using some a priori 
information e.g. [5-7] or by developing non iterative procedures. Some of these methods [6, 7] use a 
factorization approach while others are based on reformulating the inverse problems in terms of integral 
equations [8-10]. 

In [9] we presented one such approach which we applied to a two dimensional bounded domain while 
in [10] we extended this work to an unbounded three dimensional region. A feature of this method 
is the use of a priori information to find the regularized conductivity. This represents an additional 
requirement of our method but in some cases information on the internal structure of the object is 
already available. In return this enables us to make an estimate of the conductivity using a single 
data set. Another feature is that many of the calculations involve working with analytical expressions 
containing the eigenfunctions of the kernel of these equations, the computational part being restricted of 
the introduction of the data, the numerical evaluation of some of the analytic formula and the solution 
of a final integral equation. However, this first version of our method requires the knowledge on the 
boundary d^l of the region not only of the electrical potential <I> and its normal derivative d^/dn, but 
also of the electrical conductivity a and its normal derivative da/dn. In geophysics the conductivity 



function at the surface can be measured by taking small soil samples, whilst in medical applications this 
might be achieved using closely spaced electrodes. Although in geophysics it is sometimes also possible 
to measure the derivative of the conductivity, in medical applications this is more difficult. 

In this paper we describe a new version of our approach which no longer requires the knowledge of 
da/dn, but uses only data which is easily measurable on the boundary: namely the boundary values of 
the conductivity function cr, the injected currents ( i.e. d(^/dn\Q^) and the corresponding values of the 
induced potentials We will show that this information is sufficient to find a regularised solution 

Gregi^) for the Conductivity. This new approach differs from our previous work [9,10] also in that the 
final step in the reconstruction algorithm involves the solution of a first order linear partial differential 
equation rather than one of second order. 

Throughout this paper we shall assume that the conductivity distribution a G C^(i7). Such functions 
are dense in the larger class of the functions which have these properties only piecewise and which 
contains many of the physically interesting situations. In addition we shall assume that the domain J7 
is a bounded open set whose boundary dO. is sufficiently smooth, namely of class C^. 

Identification problems for elliptic equations have been the object of many studies in other fields as 
hydrology and diffusion equation [11-15]. 

2. Problem formulation 

From the conservation of the current density j(x) it follows that in the case of a general isotropic 
conductivity distribution cr(x), the potential $(x) satisfies : 

V • [o-(x)V$(x)] = 0, xGJIcM", n = 2,3. (1) 

The inverse conductivity problem considered in this paper can be stated as follows: from the values of 
(7, $ and d^/dn on the boundary dVt, reconstruct the conductivity function (^(x) in Q,. 

If we define ct(x) = ln{a{x)), equation can be rewritten as 

V2$(x) = -Y{x) , with Y(x) = Va(x) • V^'(x) . (2) 

This equation has the advantage that the {not singular) solution cr(x) to the inverse problem formulated 
in this way ensures the positivity of the conductivity (t(x). To find cr(x) we will first reconstruct the 
potential $(x) everywhere in from the information available on the boundary. We will then solve the 
differential equation for (t(x) using the method of characteristics. In this case, equation Q reads 

where s is the distance measured along the characteristics. 

By using Green's formula, we can transform the partial differential equation Q and the boundary 
information into an equivalent pair of integral equations 

cl>(x) = Xd(x) + y"d"yg^(x,y)y(y) (4) 

n 

<I>(x) = X7v(x) + J d"ya7v(x,y)y(y), (5) 

where Qo and Qn are the Dirichlet and Neumann Green's functions for Laplace's equation in Q. Here 
Xd ^-nd xn are the two harmonic functions constructed respectively from the measured Dirichlet and 
Neumann boundary data, ^\qq and d^/dn\QQ : 

Xd(x) = -J d"-iz^(x,z) Hz)\g^ (6) 
an 
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Xiv(x) = j d"-izgiv(x,z) 



5$ 
dn 



dQ 



where C = J ^{z)(P ^z. These functions are different unless o"(x) is a constant. 

dn 

Theorem 1 T/ie knowledge of cr{x) on d^l and o/<I>(x) in Q uniquely determines 0"(x) in Q. 

Proof Suppose that cti and CT2 are two functions satisfying equation © and ai / a2 • Define a^j // = 
CTi — a2 • Then 

V2$(x) = -Vai(x) • V$(x) = -V52(x) • V$(x) , x G (8) 

V^rfi//- V$(x) = 0, XGJ7. (9) 

In other words, the function a^iffi'^) is constant along the current lines inside J7. Since a is known 
on we have adiffi:^.) = 0, x G Hence, provided that through any point x G there passes a 
current line that intersects the boundary dfl [see remark below and [12]], it follows that 

adiff{^) = 0, xGf). (10) 

Thus, a is unique and it follows from the monotonicity of the logarithmic function that so too is o". □ 

Remark It is important to verify that through each interior point of 0, passes a current hnc which continues up to the 
boundary. The danger with the current hues of a general vector field is that they may spin indefinitely without reaching the 
boundary, as happens with solenoidal fields. Since in our case the vector field is of gradient type, it is irrotational and so the 
field lines cannot terminate inside the domain [16] but must continue up to the boundary. Of course, there might be some 
singular points where, for example, the gradient vanishes but as is known from Morse Theory [17], as long we are interested 
only in generic cases, these singular points are isolated. There exist extreme situations (e.g. constant fields) for which the 
current lines avoid large regions, but any generic changes, even infinitesimal, will lead to one of these (generic) Morse Fields 
in which cr can be constructed from its boundary values everywhere, except at most at some isolated points where it can be 
found using continuity. 

3. The Direct Problem 

Before considering the Inverse Problem, let us see how one can compute the potential $ if the conduc- 
tivity cr is known. For instance, as shown in [9,10], one may use the well known change of variable 
r = ^/o^ to transform equation into an integral equation for the function ^ = r<I>. Another way 
consists of applying the operator VxCr(x) • Vx to equation @ to give 

y(x) = Vx5(x) • VxXi?(x) + I d"yi^(x,y)r(y), (11) 

n 

where K{x,y) = Vxcr(x) • Vx^d(x, y). This is an integral equation for the Laplacian ^(x) of <&(x). 
Once Y{x.) is known, one can compute $(x) by means of a quadrature using the formula (0]) or ©. 

It would be nice for the practical solution of equation (|11() if the integral operator defined by the 

kernel i^r(x,y) were compact, since compact operators imply the Fredholm Alternative [18]. For that 

it is sufficient that K G L^(O^), but / d'^x J d^y |VxCr(x) • Vx^D(x,y)|^ is divergent because of the 

n n 

singularity of Vx^d(x, y). However, we shall show that the situation can be recovered by considering 
the first iteration of equation ((TT|l . 

Theorem 2 // cr G C^{Q), the first iterated kernel of equation mi]) i^2(x, z) = f d"y i^(x, y)i^(y, z) 

n 

is Hilbert- Schmidt. 



3 



Proof We shall prove that the first iterated kernel 

i^2(x,z) = y (iVVxa(x)-VxaD(x,y)Vy5(y)-Vygi5(y,z) (12) 
n 

is Hilbert-Schmidt. For smooth conductivities a in fi, any singular behaviour in the integrand comes 
from the gradient of the Dirichlet Green's function. We consider here only the three dimensional case 
but the proof in two dimensions is similar. Although in two dimensions the first iterated kernel K2 has 
a weak singularity (of logarithmic type) it is also Hilbert-Schmidt. 

We study only the leading singularity which is of the form 

Vx^D(x,y) ~ — I 13 . 

47r |x — yl'^ 

Since a E C^{Q), for any e > there exists a ball Qji C of radius R with centre x G J7 such that 

|Vx5(x) - Vy5(y)| < e , yyenR 

Hence, in this ball we can approximate Vy?(y) by Vx(t(x). The most singular parts of ^2(x, z) can be 
written as 

i^2(x, z) ~ K^{^, z) + E:2(x, z) (13) 

where 



n\nii 

and 



|y - zP 



i^|(x,z) = / dVVx5(x) . f^Vx5(x) . 

lovr^ J |x — yl'^ 



Since the modulus of the gradient of a is bounded on fi, we have 

Const /" ,0 1 



|y - zp ' 



\T^lr M Const /■ 3 
I^2(X,Z)|<^^ J d% 



X — yP |y — zP 



n\nji 



When studying the limit r = |x — z| — > we can take r < R/N « R for some N & N. Then, for 
y G it follows that |y — x| > R and |y — z| > R{1 — 1/N), and hence 

I T^i / Const 101 

LfCa(x, z) < ^- ' , ' ^ , . (14) 

To evaluate KKx, z) we define x = (0,0,0), z = (0,0, r) and y = (/5sin6'cos0, psin6'sin(/), /9cos6') , 
where {p, 9, (/)) G R'^ are local spherical coordinates in The constant gradient of a points in an 
arbitrary direction, so for illustration we shall consider Vx5(x) = (0, 0, 1). In this case 



R TT 

ryZf ^ ^ f r f ■ nrn COS 9 (o COS 9 — v) 1 . ^. 

K^ix,z) ^ — dp sm9d9 = -rpr-^- (15) 

(p2 -2prcos0 + r2)2 127ri^ 



Svr 







Combining equations (|14() and H15() we see that -fC2(x, z) remains finite as |x — z| = r — > and hence is 
Hilbert-Schmidt. □ 

This theorem implies that the first iteration of equation (|lip . namely 

y(x) = Vx?(x) • VxXd(x) + J d"y i^(x,y)Vy5(y) • VyXD(y) + J d"y ^2(x, y)) F (y) (16) 
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can be solved by means of the usual numerical procedures for Fredholm equations. Once the function 
y(x) has been computed, the potential ^*(x) can be obtained by means of the formulae @ or ©. This 
procedure can then be used to perform simulations to study the performance of our method of solving 
the Inverse Problem. 

4. The Unregularized Inverse Problem 

Initially, our treatment of the inverse problem is similar to that described in [9,10]. Subtracting equation 
(jlj) from © we have 

X(x)- J d"y/C(x,y)y(y)=0 (17) 

n 

where 

^(x,y) = ^D(x,y) - ^Ar(x,y) and x(x) = Xn(x) - xd(x) . 

Eauation H17(l is an ill-posed Fredholm integral equation of the first kind. Since the kernel /C(x,y) is the 
difference of two Green's functions it is harmonic and symmetric. The null space N(IC) of /C consists 
of all the functions which are orthogonal to the harmonic functions in Q. We can prove the following 
result concerning the potentials $j related to this null space. 



Theorem 3 Within the set of functions {^^} that satisfy V^<I>j = —1^ where G ^{^) there exists 

0. 



at least one function such that <I>j(x)|g^ = and -^(x 



an 



Note: Such a function is, of course, completely invisible to any measurement of potentials and currents 
on the boundary. 

Proof Since the Laplacian of any harmonic function is identically zero, it is clear that the functions 
$j are defined only up to an harmonic function ^h- So, if the boundary values of some initial <I>^ are 
not zero, we can always find another function <I>j = <I>^ + having the same Laplacian (— Yj), but 
with 

^t\an = 0- (18) 

To this end it is sufficient to take 

4..(x)=/d"-'y^?£^*i(y). (19) 

act 

where the function $/j is harmonic and has by construction boundary values which are opposite to that 

of^r So ^i\s^ = o. 

To prove that the normal derivative of $j is identically zero on the boundary, we use Green's formula 



j ('u(x)V2«>j(x) - $5(x)V2n(x)) = j (f"-ix (^'u(x)^^^ - «>j(x) 



an 

where u(x) is any harmonic function in Q. Since V^u = and since V^<I>| = —Y^ is orthogonal to any 
harmonic function, the left hand side of the above equation is zero. But from H18() 'I'jIgQ = 0, we see 



that 

dn 



I 



<f-'x„(x)52jM.„^ (21) 



an 

Let us assume for a moment that at some given point xq G the normal derivative were not zero, 
5$j(x)/9n|^^^^ 7^ 0. Since according to our general assumptions the derivatives of $ are continuous 
in — i.e. also on 50 — there exists a neighbourhood on the boundary dVt^xo C dVt of xq S 90, 
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on which this derivative has the same sign. Considering for u(x) the harmonic measm'e defined by the 
boundary conditions 

u[x)\^(,dn - I , elsewhere on dQ. ^ > 

we see that if d^^/dn were not identically zero on dQ, the integral H21() cannot be zero. As a consequence 
both $j(x) and d^^{x)/dn have to vanish identically on d^l. □ 

An example of such a function can be constructed explicitly from the function Y^, as described in the 
following Lemma. 



Lemma 1 The function ^'j(x) = J (i^y ^£)(x, y)Yj(y), where £ N{JC), has the properties stated 

n 

Theorem 3 i.e.: (i) V'^<^^ = -Y^ x e n, {ii) ^>j(x) = , x E 90, (Hi) ^(x) = , x e dn . 
Proof (i) If X G r2 , then 

V2$e(x) = j d"yV2gB(x,y)yj(y) = - j d"y5(x - y)yj(y) = -yj(y) . 
n n 

(a) If X G dn, recall that C/D(x,y) = and hence ^^ix) = J (i"y ^£i(x,y)yj(y) = 0. 
(Hi) If y^ G N{IC) and /C(x, y) = ^d(x, y) — ^Ar(x, y) it follows that 

$e(x)= I d"yaiv(x,y)yt(y), xgO 
is also a valid representation for $6. Moreover, ^^^(^'y) 



— I 1 so that 



9$j 

9n 



(x) 



\dn\ 



since (y) is orthogonal to any harmonic function in 0, and the constant 1 is an harmonic function. □ 

Explicit analytical expressions of such functions Y^ G N{IC) and their respective <I>j for the unit disk are 
given in Section 6. 

At this stage an apparent paradox may emerge. Consider the case $ = corresponding to Y^ G N{IC). 
Since is in general not identically zero throughout 0, (see, for example, the explicit expressions 
given in Section 6), one may be tempted to reconstruct a inside 0, by solving along the characteristics 
the differential equation ^ using the given boundary data (t(x) = \n{a{x)) for x G dn. However, 

J (T'xaix) |V«>(x)|2 = j d"-ixa(x)$j(x)^^^ =0, (23) 



d<S>r, 



since both $^( x)|qq and (x) are identically zero. From the non-negativity of (t(x) — exp(cj(x)), 
the left hand side is positively defined and so a has to be zero everywhere where V<I>j is not. The key of 
this apparent dilema is that equation Q cannot be solved since the normal component of V<^j vanishes 

when X approaches the boundary. This represents a barrier that cannot be overpassed ( ^ oo) by 
the information conveyed by the equation Q. On the other hand, if a is constant, from equation ((2)) it 
follows that y^ has to vanish identically throughout n, i.e. the coefficients A^^n and Bk^n in expression 
have to be zero. 
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5. The Regularized Integral Equation 

Equation (|17jl is a Fredholm equation of the first kind which is notoriously ill-posed and whose solution is 
extremely unstable with respect to the data measured on the boundary Ofi. We shall stabilize it by using 
Tikhonov regularization and a model function Ymodi^) as a first estimate of the solution (see Figure 

This function can be obtained starting with a model conductivity amod and solving equation 
or better its first iterate (fT6|) . to find Ymod(x). The function amod can be obtained in a variety of ways 
including the use of prior knowledge from other imaging modalities, previous results in time evolving 
situations or using simple reconstruction estimates using fast algorithms. One example of the latter is a 
Factorization method described in [6,7]. In this method a large number of current patterns are applied 
to the object and using this data a simple test at each point within O is used to determine a piecewise 
constant internal structure. This algorithm gives good results in many circumstances although practical 
numerical difficulties introduce instability. Nevertheless, we can use a continuous approximation to these 
results as a suitable amod- 

The use of the regularizing function Ymod leads to a Fredholm equation of the second kind which can 
be solved explicitly and which is stable with respect to the errors in the input data. Moreover the null 
space of /C plays no particular role in the solution of Fredholm equation of the second kind which is 
unique unless the Lagrange multiplier defined below happens to be an eigenvalue of IC. 

Indeed, the input data and hence the function x(x) always contain some errors. Consequently we no 
longer require that equation ()17|) should be satisfied exactly, but we ask instead that the norm of 
the left hand side of (fT7|) should be small while the solution Yreg should be close to a model function 
Ymod- Hence our problem reduces to finding 

eo = min\\x{x) - K[Yreg]{:>i) Wl^ , subject to ||l"reg(x) - Ymod{x)\\L2 < S. (24) 

Here K[yreg](x) = / d^y IC{x,y)Yregiy)- In practical terms this means that instead of eq.(|17|) we 
consider the Lagrange multiplier formulation of Tikhonov regularization [19,20], namely we seek the 
unrestricted minimum of the functional 

IIX(X) -K[Yreg]{^)\\l, + /i (^.^^(x) - ^^(x)!!^^ " S^) , 

where is a Lagrange multiplier. 

In many cases of practical interest the target function depends also on a small number of parameters a, 
for instance connected with the translation and/or rotation Euclidean group. Consequently we could 
include these parameters in the model function l^od(x; a) and seek the minimum eo^min = eo(cKo) with 
respect to these parameters. Such an example is discussed in Section 8. Here eo which is a functional 
of the data, or its minimum eo,mjn with respect to the possible parameters a, measures the degree [20] 
to which the initial integral equation (|17|) is satisfied. In other words, it reflects the accuracy of the 
experimental data by means of which (see eq.lO) the input function x(x) has been constructed. 

The minimisation process yields the following regularized integral equation: 

Yregi^) = Ymod{^) + X J d"y /C(x, y) x(y) - \ J (i"y ^2 (x, y) F^e^ (y) (25) 

n n 

where /C2 (x, y) = J d"z /C(x, z) /C(z, y) andA = l//i. 

n 

The eigenfunctions {ufc} and eigenvalues {A^} of /C(x, y) satisfy Ufc(x) = Xk f (^"y ^(x, y) itfc(y). If the 

n 

kernel /C(x, y) can be expressed in the form 



^(x,y) = 2^- 



U 1 

k=l 



7 




Figure 1: 

For the Tikhonov regularisation it is not necessary that the mapping A4 between the parameter space C 

and the data space S should have a true inverse. The existence of a generalized inverse is sufficient. 

In the above figure the mapping transforming the C plane into the straight line Ac G <S is provided by the 

non-invertible matrix ^ ^ while the generalized inverse transforms the points Pi, Po G Ac into 

the lines Di, Do G C. 

These low dimensional graphs may be misleading, as in infinite dimensions it might happen that there does 
not exist a smallest distance eo between the admissible set Ac = MiC) and the data ^ Ac- Therefore, in 
infinite dimensional problems it is preferable to consider images of balls B from C rather than M{C) itself. 
Since in most physical problems the theorem of Banach-Alaoglu [22,23] also applies [24] , the images A4{B) 
are true compact sets and the minimal distance eo of such a set to the data is well defined. Nevertheless, 
the main features of these problems are well depicted by these figures, where the interesting points on the 
straight lines Di ...,Do are their intersections with the ball Bs of radius S around Ymod- Consequently, in 
our case the solution of the extremal problem will be the point Yreg G Do. 

we have a straightforward way to compute the eigenfunctions and eigenvalues and we can expand the 
iterated kernel IC2 as: 

J. . . ^ Ufc(x) Ufc(y) 
k=i 

This is illustrated in the two dimensional example discussed in Section 6 (see also [21]). 
Projecting Eq.((^SJ) onto Uk we obtain 

^reg,k ^mod,k ~l~ T ~X.k T2" ^reg,k 



where 



Yreg,k = j (f'y^Yreg{^)uk{y^) , Ymod,k = j (i"x ymod(x)lifc (x) , Xk = j (i"x ^(x)^^; (x) , 



n n 
and hence the following explicit expression for the solution Yreg{x.): 



Ymod,k + T— XA: 



dy" /C(x, y) x(y) - A ^ ^ ^— u^i^) . (26) 

II A + Ai. 



Once l^eg(x) is known, the potential $reg(x) can be obtained by means of the integral 

$reg(x) = Xd(x) + J d"y (x, y) F.e^ (Y ) • (27) 

Finally, to determine the regularized (Treg(x), we can use the method of characteristics to solve the first 
order partial differential equation: 

Varegix) ■ V^regi^) - Yreg{^) = , X G (28) 
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subject to the known boundary values (^{^)\xedn- 

6. A two dimensional example: the unit disk 

In this case the relevant Green's functions are: 

^ , ^ „N 1 , + p2 _ 2rpcos (6 - 1}) 

GDir, e;p,^) = - — log / ^ 29 

47r 1 + — 2rp cos [a — v) 

gN{r,e;p,'d) = --^ logUr'^ + p'^ -2rp cos {6 -■&))■ p"^ -2rp cos {6 -'&))) (30) 

47r 

so that 

X:(r,0;/,,z?) = — log (l + rV-2rpcos (0-19)) - ^V" cos(fc(g - ^)) _ ^^^^ 

27r 1TK 

k=l 

As remarked in Section 5, this series expansion enables us to find the eigenfunctions and eigenvalues 
of /C. Indeed, the functions u^(r, 6) = {2k + 2)/7r cos kO and u\{r, 9) = \/ {2k + 2)/7r r*" sin /c^ are 
orthonormal on the unit disk, so that /C can be rewritten as 

k=l j=l ^ 

with Afc = -2fc(A; + 1). 

From expansion 1)26^ we obtain the following explicit expression for the solution of the regularized 
equation 

„ yj ^ 

/■I /"Stt °° ^ mod,k OA'i'fr -I- 1 

Yre,{r,e) = Y^Ur,e) + \ j pdp j IC{r,9; p,^) x{p,^) - X Yl E a + 4P(fc + l)2 

(32) 



To find dreg we proceed in the way discussed at the end of Section 5. 
Lemma 2 For the unit disk the functions G ^{^) have the explicit form 



oo oo 

^«(^' ^) = E E '^"(^ + 2, A; + 2, r) (a^,, cos(A;e) + S^,, sin(A:e)) , (33) 



fc=i ji=i 

where Aj,^n c-i^d B^^n (^^^ arbitrary constants. 
Proof If e iV(/C), then 



»1 »27r ^ 

/ pdp / d^/C(r, 0; p, i?)yj(p, ^) = Y.Y. 
■^0 -^0 fc=i ,=1 



u{{r,9^ '■2- 



We start with a function of the form 

Y^{r,6) = ^k{r)cos{k9), k = 1,2,3, 
which, in order to be in N{IC) must satisfy 

2tv 



pdp d^ui{p,^)Yi{p,'d) = 0. 
Jo 



Jo 



pdp / d^ul{r,9;p,^)Y^{p,^) = -^{2k + 2)7r / dp p^+^^k{p) = 0. 
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For this it is sufficient that 



3f?fc(r, n) = Gn{k + 2, + 2, r) , forn > 1 , 

where Gn{k + 2, A; + 2, r) are the Jacobi polynomials of order n [18]. 
If we consider 

yj(r, 9) = 3?fc(r) sin{ke) , A; = 1, 2, 3, ... . 
we find a similar result. Our theorem follows by taking linear combinations of these results. □ 
Example: If we consider Yj(r, 9) = d (3, 3, r) cos(6l) = (l - |?^) cos(6l) G A^(/C) we can compute directly 

^l{r,e)= j\dp j^J d^GD{r,9-p,'&)Y^{p^) = h{r - if cos{9) . 



7. A numerical example 

We present here a numerical example to illustrate the performance of our algorithm. In this example 
we attempt to reconstruct a conductivity distribution consisting of two low conductivity regions and a 
high conductivity region within an uniform background. We represent this distribution mathematically 
by a function of the form: 

cTexact(x, y) = I - 0.5 exp(-a((a; - 0.5)^ + (y - 0.2)2)2) ^ exp(-6((x + 0.1)^ + (y - 0.3)2)2) 

- 0.5 exp(-a((x + 0.5)2 + + 0.2)2)2) ^ + y'^ <l, (34) 

which has a maximum at {x,y) = (—0.1,0.3) and two minima at(x,?/) = (—0.5,-0.2) and (x,y) = 
(0.5,0.2) (see Figure |2(^. We chose a = 1000 and h = 1500. 



We shall consider the input currents 



a{l,e) 



dn 



{9) = sm(m9) , m = 1, 2, . . . . 



an 



To simulate the measured values of the potential on the boundary we have to solve first the direct 
problem. The detailed numerical simulations can be found in [25]. With these data we can compute 
then the input function x of our integral equation 



Xir,9) 







2^^ r d<^ 

d^lgN{r,9;l,^) — 



an 



27r 



dn 



^>(1 



To calculate the regularized solution Y^eg we need a model function Y^od- To achieve this we need 
a model function for the conductivity amod^ which is in fact our a priori information. For the model 
function (Jmod we have used a conductivity distribution in which 

- the positions of the low conductivity regions are quite well known but the size of these regions is more 
uncertain 

- the size of the high conductivity region is quite well known but the position of this region is uncertain. 
We have represented this model function mathematically as: 



(^mod{x,y) 



1 - 0.7exp(-a((x - 0.6)2 + (y _ 0.3)2)^) + 1.2exp(-5((x + 0.4)2 + _ 0.7) 



\2\2\ 



0.7exp(-a((a; + 0.6)^ + (y + 0.3)2)^) , 



(35) 



which lies at a distance d = WcTmod ~ ^exactW l'^ 

= 0.2911 

Inexact 1 1 l2 from the true function (see Figure [2 (b)| ). 
This type of function could be a smooth approximation to piecewise constant approximates produced 
by [6, 7] as discussed in Section 5. By solving the forward problem as described in Section 3, we find 
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(a) (b) 
Fi gure 2: (a) The exact conductivity Cexactj (b) The model conductivity Cmod- 

the model potential ^mod and hence the model functions Ymodix,y) = Vamod{x,y) • V^mod{x,y)- By 
solving the forward problem as described in Section 3, we find the model potential ^mod and hence the 
model functions Ymodix,y) = Vo'modix,y) • V^mod{x,y). The integral which appears in (|32|) can be 
computed numerically: 

/ pdp / di?/C(r, 6; p, ^)xip, ^) = V WifCir, 9; n,0i)xin,ei) , 
Jo Jo 

where {wi;ri,6i} is a set of quadrature weights and points for the unit disk [26]. Since, the infinite 
series which appear in equation (|32|) converges rapidly it was sufficient to consider only fifty terms to 
reach a precision of 10~®. 

There are many ways of determining the value of the parameter A [19,27] and we have used the constraint 
(|24jl for a given value of 5. Once Yreg is known, the calculation of the potential $ inside the unit disk 
is straightforward: 

172 

^{r, 6) = XD{r, 0) + Y, ^lOoir, 9; ri,ei)Yreg{ri,ei) , 
1=1 

where 

^^^"'^^-^^X l + r2-2rcos(e-^)^^- 

The final step is to solve our first order partial differential equation (|28|) . The method of characteristics 
transforms this equation into an equivalent system of ordinary differential equations which can be 
written in cartesian coordinates as: 

dx d^ dy d^ da 

- = -(x(.),y(.)), ^ = ^(^(^),y(«)), ^ = YreMs)Ms))- (36) 
We have used sixty characteristic curves, the initial conditions for the i-th characteristic being: 

Xo = cos0O) yo = sin^oi ^o = 27r/i, /or i = l... 60, (37) 

5^ = ln(l - 0.5 exp(-1000(1.29 - cos O'^ - 0.4 sin e'^f) + exp(-1500(l.l + 0.2 cos 6*^ - 0.6 sin 6*^)2) 
- 0.5exp(-1000(1.29 + cos (9^ + 0.4 sin (9o)^)) . 
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We present in Figures |31 and 0] the results obtained for two current patterns with m = 1 and 2 for data 
with 2% random errors. The reconstructed conductivity obtained hes at a distance \\a — (7exact\\L^ = 
0.07||(Jexarf|lL2 for 171 = 1 and at \\a — (TexactWh"^ = 0.09||(Tei:act||L2 for m = 2 from the true function. We 
see that in the both cases the norm of the difference between the target and the model has been 
reduced significantly. It is interesting to note that these results are quite similar and that it appears 
that little can be gained by using more patterns. In our method the interest of using extra current 
patterns is in the calculation of a good model function Ymod- 





8. The notch of the eo graph 

We finally discuss an example, in which the position and/or orientation of the target function is described 
by a small number of parameters. Such a situation occurs for example in land mine detection, where 
the type of objects we are seeking (i.e. the model function) is quite well known and it is their location 
which we wish to determine. This example takes into consideration the action of some group of operators 
acting on the model function. In practice we shall encounter mainly the action of the Euclidean group 
of rotation and translations. 
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An interesting feature of this kind of problem is that we are no longer in 
the ball regularization geometry where one often encounters the situation 
in which the final solution does not represent an improvement with respect 
to the model function. This arises in simple ball regularization since the 
fact that both the target function Yexact and the solution l^eg lie on the 
surface of an infinite dimensional sphere of radius S centered at Y^od does 
not necessarly imply that \\Yreg - YexactWi^ < \\Ymod - YexactWh^ = 5 (see 
Figure • In two dimensions we can see that the probability that this is 
true is only 1/3 as Yexact needs to lie on the 120" bold face arc around Yreg, 



Ymod 



Figure 



Yreg 

5: In two di- 

while, in the infinite dimensional case this probability reduces to 1/4. This mensions the probability 

that litres - YexactWL'i < 

\\Y„iod - YexactWL^ = 5 IS Only 
1/3 as Yexact needs to lie on the 



is certainly disappointing, but this situation will no longer hold if Ymod 
depends on some unknown parameter(s) a so that the geometry of the 
regularization set is no longer spherical. We then try to find the optimal 
value ao of the parameter a (see Section 5) by plotting the graph of eo(a) 
defined in equation H24() . 



120° bold face arc around Yn 





(a) The translation of amod towards (Jexact for different values 
a. 



(b) The notch of the graph eo («) for 
two different 5 = \\Yreg - YexactWL^ 



Figure 6: The notch method. 



As an example, we consider the problem of the lower half plane in which the conductivity depends on 
a parameter a. The target conductivity has a maximum on the y = —1 axis at a point corresponding 
to ao = 1.5. To illustrate how this procedure works we consider only the translations of the model 
function along the line y = —1, for a G [0,3] (see FigureEl (a)). We represent the target conductivity 
distribution by 

— ^ 

~ ^ ((x-ao)2 + (?/ + l)2)2 + 6 

and we choose the model function of the form 

CTmodei = 1 + exp(-c((a; - af + {y + iff) ■ 
We chose a = b = 0.001 and c = 1500. 

By means of a conformal mapping we transfer the problem onto the unit disk and proceed as described 
earlier for different values of a. Since this is only an illustrating example, in the regularization functional 
(|24)1 we used the unweighted norm for the unit disk which means that for the initial half plane we 
used the norm induced by the Jacobian of the inverse mapping. This norm has compactifying properties 
but it depends on the conformal mapping. Even in this rough treatement we see that this method is 
quite sensitive since the notch in the graph of eo(Q^) is quite narrow (see Figure |H1 (b)). The position of 
this minimum of the graph eo(a) gives the location of the object we look for and represents the best 
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(a) The regularized conductivity are.g- 



(b) The absolute errors Icrreg — creaiactl- 



Figure 7: The reconstructed conductivity using the notch method. 

value of a to use in our model. The conductivity distribution reconstructed using this method is shown 
in Figure [71 

9. Conclusion 

We have presented a reconstruction algorithm based on a Tikhonov regularization of a formulation of 
the inverse conductivity problem in terms of integral equations. The method uses prior information in 
the form of a model conductivity distribution. Our numerical experiments have shown the importance 
of this model and that multiple data appear to have little effect on our results. We have suggested ways 
in which more data measurements may be used to improve the model but numerical results suggest that 
our technique may prove to be useful in situations in which complex data collection is difficult but we 
have a good qualitative information on the conductivity distribution . 
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